import numpy as np
import scipy as scp
import matplotlib.pyplot as plt
plt.style.use(['dark_background'])
import sympy as sp
import sympy.physics.hydrogen as sphydrogen
import scipy.constants as const
$E_h = 27.21138602 \,eV\, = 4.35974434 × 10^{-18} \, J$
I have defined the functions in the Defining Function section. Here I have taken functions from sympy.physics.hydrogen and created numerical functions from them. The defined function, named hydrogen_wavefn_all gives the analytical and graphical view of total wavefunction $\psi_{nlm}(r, \theta, \phi)$ and radial wavefunction $R_{nl}(r)$ and hydrogen_wavefn_return returns 2 functions phi_nlm_fn(r, theta, phi) and R_nl_fn(r) for those wavefunctions.
In the Give Inputs and Get States section, you can give arbitary values of (n, l, m) and get the states. Here you need to set the number of grid points (N1) and the maximum limit of radius r1 to view the total graphs.
Once you are done with the above things, go to the Working with Functions section. Here I have used the returned functions (phi_nlm_fn(r, theta, phi) and R_nl_fn(r)) and done few more plots from these.
hydrogen_wavefn_all¶cmap: 'Accent', 'Accent_r', 'Blues', 'Blues_r', 'BrBG', 'BrBG_r', 'BuGn', 'BuGn_r', 'BuPu', 'BuPu_r', 'CMRmap', 'CMRmap_r', 'Dark2', 'Dark2_r', 'GnBu', 'GnBu_r', 'Greens', 'Greens_r', 'Greys', 'Greys_r', 'OrRd', 'OrRd_r', 'Oranges', 'Oranges_r', 'PRGn', 'PRGn_r', 'Paired', 'Paired_r', 'Pastel1', 'Pastel1_r', 'Pastel2', 'Pastel2_r', 'PiYG', 'PiYG_r', 'PuBu', 'PuBuGn', 'PuBuGn_r', 'PuBu_r', 'PuOr', 'PuOr_r', 'PuRd', 'PuRd_r', 'Purples', 'Purples_r', 'RdBu', 'RdBu_r', 'RdGy', 'RdGy_r', 'RdPu', 'RdPu_r', 'RdYlBu', 'RdYlBu_r', 'RdYlGn', 'RdYlGn_r', 'Reds', 'Reds_r', 'Set1', 'Set1_r', 'Set2', 'Set2_r', 'Set3', 'Set3_r', 'Spectral', 'Spectral_r', 'Wistia', 'Wistia_r', 'YlGn', 'YlGnBu', 'YlGnBu_r', 'YlGn_r', 'YlOrBr', 'YlOrBr_r', 'YlOrRd', 'YlOrRd_r', 'afmhot', 'afmhot_r', 'autumn', 'autumn_r', 'binary', 'binary_r', 'bone', 'bone_r', 'brg', 'brg_r', 'bwr', 'bwr_r', 'cividis', 'cividis_r', 'cool', 'cool_r', 'coolwarm', 'coolwarm_r', 'copper', 'copper_r', 'cubehelix', 'cubehelix_r', 'flag', 'flag_r', 'gist_earth', 'gist_earth_r', 'gist_gray', 'gist_gray_r', 'gist_heat', 'gist_heat_r', 'gist_ncar', 'gist_ncar_r', 'gist_rainbow', 'gist_rainbow_r', 'gist_stern', 'gist_stern_r', 'gist_yarg', 'gist_yarg_r', 'gnuplot', 'gnuplot2', 'gnuplot2_r', 'gnuplot_r', 'gray', 'gray_r', 'hot', 'hot_r', 'hsv', 'hsv_r', 'inferno', 'inferno_r', 'jet', 'jet_r', 'magma', 'magma_r', 'nipy_spectral', 'nipy_spectral_r', 'ocean', 'ocean_r', 'pink', 'pink_r', 'plasma', 'plasma_r', 'prism', 'prism_r', 'rainbow', 'rainbow_r', 'seismic', 'seismic_r', 'spring', 'spring_r', 'summer', 'summer_r', 'tab10', 'tab10_r', 'tab20', 'tab20_r', 'tab20b', 'tab20b_r', 'tab20c', 'tab20c_r', 'terrain', 'terrain_r', 'turbo', 'turbo_r', 'twilight', 'twilight_r', 'twilight_shifted', 'twilight_shifted_r', 'viridis', 'viridis_r', 'winter', 'winter_r'
Selected: 'viridis', 'plasma', 'inferno', 'magma'
def hydrogen_wavefn_all(Z, n, l, m, r, theta, phi):
'''
------- Functions to import --------
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
import sympy.physics.hydrogen as sphydrogen
-------- Inputs ---------
Z: atomic number
n: principal quantum number: 1,2,3,...
l: azimuthal quantum number: 0,1,2,... (l<n-1)
m: magnetic quantum number: ...,-2,-1,0,1,2,... (|m|<l)
r: an numpy array: np.linspace(0,15,N)
theta: numpy array: np.linspace(0, np.pi, N)
phi: numpy array: np.linspace(0, 2*np.pi, N1)
N = 50 # number of grid points for plotting
------- Outputs ----------
value of energy of the state
Analytical expression for psi_nlm(r,theta,phi)
Normalization check of psi_nlm(r, theta, phi) (compare with 1)
Contour plot for psi in (r*cos(psi), r*sin(phi)) plane
Analytical expression for R_nl(r)
Normalization check of R_nl(r) (compare with 1)
Plot of r vs R_nl(r)
'''
print(f'Inputs: n={n}, l={l}, m={m}, Z={Z}')
rsp = sp.Symbol('r', positive=True)
phisp, thetasp = sp.symbols('phi theta', real=True)
Zsp = sp.Symbol('Z', positive=True, integer=True, nonzero=True)
psinlmsp, rnlsp = sp.symbols('\psi_{nlm} R_{nl}')
EheV, EhJ = 27.21138602, 4.35974434e-18 # Energy in eV and J
E_n = sphydrogen.E_nl(n, Z)
print(f'E = {E_n} in Hartree atomic units = {EhJ*E_n} J = {EheV*E_n} eV')
Psi_nlm_sp = sphydrogen.Psi_nlm(n, l, m, rsp, phisp, thetasp, Zsp)
display(psinlmsp, Psi_nlm_sp, Psi_nlm_sp.simplify())
print('normalization:', sp.integrate(sp.Abs(Psi_nlm_sp)**2*rsp**2*sp.sin(thetasp),
(rsp,0,sp.oo),(thetasp,0,sp.pi),(phisp,0,2*sp.pi)))
Psi_nlm_np = sp.lambdify((rsp,phisp,thetasp), sphydrogen.Psi_nlm(n,l,m,rsp,phisp,thetasp,Z))
# R3, Theta3, Phi3 = np.meshgrid(r, theta, phi, indexing='ij')
# psiarr3 = Psi_nlm_np(R3, Theta3, Phi3)
R2, Phi2 = np.meshgrid(r, phi)
Theta2 = np.full_like(R2, np.pi/2) # input value of theta
psiarr2 = Psi_nlm_np(R2, Theta2, Phi2)
plt.contourf(R2*np.cos(Phi2), R2*np.sin(Phi2), psiarr2, cmap='inferno')
plt.colorbar(label='$\psi_{nlm}$')
plt.title(f'State {n, l, m}, Z={Z}')
plt.show()
plt.contourf(R2*np.cos(Phi2), R2*np.sin(Phi2), np.abs(psiarr2)**2, cmap='magma')
plt.colorbar(label='$|\psi_{nlm}|^2$')
plt.title(f'State {n, l, m}, Z={Z}')
plt.show()
R_nl_sp = sphydrogen.R_nl(n, l, rsp, Zsp)
display(rnlsp, R_nl_sp, R_nl_sp.simplify())
R_nl_np = sp.lambdify(rsp, sphydrogen.R_nl(n, l, rsp, Z))
R_nl_plot = R_nl_np(r)/(np.sum(R_nl_np(r)**2*r**2*(r[1]-r[0])))**0.5
print(f'normalization: {np.sum(R_nl_plot**2 *r**2 *(r[1]-r[0]))}')
plt.plot(r, R_nl_plot)
plt.axhline(0, color='yellow')
plt.xlabel('$r$')
plt.ylabel('$R_{nl}$', fontsize=12)
plt.title(f'State {n, l}, Z={Z}')
plt.grid()
plt.show()
plt.plot(r, R_nl_plot**2)
plt.axhline(0, color='yellow')
plt.xlabel('$r$')
plt.ylabel('$|R_{nl}|^2$', fontsize=12)
plt.title(f'State {n, l}, Z={Z}')
plt.grid()
plt.show()
print(hydrogen_wavefn_all.__doc__)
------- Functions to import --------
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
import sympy.physics.hydrogen as sphydrogen
-------- Inputs ---------
Z: atomic number
n: principal quantum number: 1,2,3,...
l: azimuthal quantum number: 0,1,2,... (l<n-1)
m: magnetic quantum number: ...,-2,-1,0,1,2,... (|m|<l)
r: an numpy array: np.linspace(0,15,N)
theta: numpy array: np.linspace(0, np.pi, N)
phi: numpy array: np.linspace(0, 2*np.pi, N1)
N = 50 # number of grid points for plotting
------- Outputs ----------
value of energy of the state
Analytical expression for psi_nlm(r,theta,phi)
Normalization check of psi_nlm(r, theta, phi) (compare with 1)
Contour plot for psi in (r*cos(psi), r*sin(phi)) plane
Analytical expression for R_nl(r)
Normalization check of R_nl(r) (compare with 1)
Plot of r vs R_nl(r)
hydrogen_wavefn_return¶def hydrogen_wavefn_return(Z, n, l, m):
'''
------- Functions to import --------
import numpy as np
import sympy as sp
import sympy.physics.hydrogen as sphydrogen
-------- Inputs ---------
Z: atomic number
n: principal quantum number: 1,2,3,...
l: azimuthal quantum number: 0,1,2,... (l<n-1)
m: magnetic quantum number: ...,-2,-1,0,1,2,... (|m|<l)
------- Returns ----------
np.array([E, E in J, E in eV]) psi_nlm(r,theta,phi), R_nl(r)
'''
# print(f'Inputs: n={n}, l={l}, m={m}, Z={Z}')
rsp = sp.Symbol('r', positive=True)
phisp, thetasp = sp.symbols('phi theta', real=True)
EheV, EhJ = 27.21138602, 4.35974434e-18 # Energy in eV and J
E_n = sphydrogen.E_nl(n, Z)
E_n_J, E_n_eV = EhJ*E_n, EheV*E_n
Psi_nlm_np = sp.lambdify((rsp,phisp,thetasp), sphydrogen.Psi_nlm(n,l,m,rsp,phisp,thetasp,Z))
R_nl_np = sp.lambdify(rsp, sphydrogen.R_nl(n, l, rsp, Z))
return np.array([E_n, E_n_J, E_n_eV]), Psi_nlm_np, R_nl_np
print(hydrogen_wavefn_return.__doc__)
------- Functions to import --------
import numpy as np
import sympy as sp
import sympy.physics.hydrogen as sphydrogen
-------- Inputs ---------
Z: atomic number
n: principal quantum number: 1,2,3,...
l: azimuthal quantum number: 0,1,2,... (l<n-1)
m: magnetic quantum number: ...,-2,-1,0,1,2,... (|m|<l)
------- Returns ----------
np.array([E, E in J, E in eV]) psi_nlm(r,theta,phi), R_nl(r)
N1 = 500
rmax = 50
Z1 = 1
n1, l1, m1 = 5,3,-1
r1 = np.linspace(0, rmax, N1)
theta1 = np.linspace(0, np.pi, N1)
phi1 = np.linspace(0, 2*np.pi, N1)
hydrogen_wavefn_all(Z1, n1,l1,m1,r1,theta1,phi1)
Inputs: n=5, l=3, m=-1, Z=1 E = -1/50 in Hartree atomic units = -8.71948868000000E-20 J = -0.544227720400000 eV
normalization: 1
C:\ProgramData\Anaconda3\lib\site-packages\numpy\ma\core.py:2829: ComplexWarning: Casting complex values to real discards the imaginary part _data = np.array(data, dtype=dtype, copy=copy,
normalization: 1.0
Main Function: hydrogen_wavefn_return(Z, n,l,m)
Returns:
np.array([E, E in J, E in eV])phi_nlm_fn(r, theta, phi)R_nl_fn(r)Z1, n1, l1, m1 = 1, 2, 1, -1 # INPUT
En, psi_nlm_fn, R_nl_fn = hydrogen_wavefn_return(Z1, n1,l1,m1)
display(En, psi_nlm_fn, R_nl_fn)
array([-1/8, -5.44968042500000e-19, -3.40142325250000], dtype=object)
<function _lambdifygenerated(r, phi, theta)>
<function _lambdifygenerated(r)>
plotly¶import plotly.graph_objects as go
N3 = 200
r3 = np.linspace(0, n1*12, N3) # input
phi3 = np.linspace(0, 2*np.pi, N3)
R2, Phi2 = np.meshgrid(r3, phi3)
theta_const = 0.5 # input
Theta2 = np.full_like(R2, theta_const)
psinlmplotly = psi_nlm_fn(R2, Theta2, Phi2)
fig1 = go.Figure(data=[go.Surface(z=np.abs(psinlmplotly)**2,
x=R2*np.cos(Phi2), y=R2*np.sin(Phi2), colorscale='turbo')])
fig1.update_layout(scene=dict(
xaxis_title='r cos(phi)',
yaxis_title='r sin(phi)',
zaxis_title='|psi_nlm|^2'),
title=f'State {n1,l1,m1}')
fig1.show()
N4 = 20
r4 = np.linspace(0, n1*12, N4) # input
theta4 = np.linspace(0, np.pi, N4)
phi4 = np.linspace(0, 2*np.pi, N4)
R3, Theta3, Phi3 = np.meshgrid(r4, theta4, phi4)
psi2nlmplotly = np.abs(psi_nlm_fn(R3, Theta3, Phi3))**2
X3, Y3 = R3*np.sin(Theta3)*np.cos(Phi3), R3*np.sin(Theta3)*np.sin(Phi3)
Z3 = R3*np.cos(Theta3)
isovalue = np.abs(psi2nlmplotly).max() *0.5
fig2 = go.Figure(data=go.Isosurface(
x = X3.flatten(),
y = Y3.flatten(),
z = Z3.flatten(),
value = psi2nlmplotly.flatten(),
isomin=0,
isomax=isovalue,
opacity = 0.7, # CHANGE
surface_count=2,
colorscale='plasma' ))
fig2.update_layout(
title=f'State {n1,l1,m1}, Energy={En[2]:.4f}eV',
scene=dict(xaxis_title='x', yaxis_title='y', zaxis_title='z'))
fig2.show()
display(psi2nlmplotly.flatten().shape, Z3.flatten().shape)
(8000,)
(8000,)
Z2, n2 = 1, 5
rmax2, N2 = (n2+2)*10, 1000
r2 = np.linspace(0, rmax2, N2)
Rnl_all = list(np.zeros(n2))
for l2 in range(n2):
En2, psinlm2, Rnl_all[l2] = hydrogen_wavefn_return(Z2,n2,l2,0)
Rnl_plot = np.array(Rnl_all[l2](r2))
print(f'normalization: {np.sum(Rnl_plot**2*r2**2*(r2[1]-r2[0]))}')
plt.plot(r2, Rnl_plot)
plt.axhline(0, color='yellow')
plt.xlabel('$r$')
plt.ylabel('$R_{nl}$', fontsize=12)
plt.title(f'State {n2, l2}, Z={Z2}')
plt.grid()
plt.show()
normalization: 0.9978348220207655
normalization: 0.9983043342534061
normalization: 0.9989981578914803
normalization: 0.9995956114985672
normalization: 0.999915281662709